function [X, L, X_err, L_err, iter, flag] = powermethod(A, Xo, X_ref, L_ref, maxit, tol, print)
%
%  [X,L,X_err,L_err,iter,flag] = power(A, Xo, X_ref, L_ref, maxit, tol)
%
%  POWER applies the power method to the eigenvalue problem
%     
%    (A - L*I)*X = 0
%
%  The iterates take the form
%    X[0]   = normalized guess, i.e. ||X[0]||_2 = 1.0
%    X[l+1] = A * X[l]
%    X[l+1] = X[l+1] / ||X[l+1]||_2
%    L[l+1] = X[l+1]' * A * X[l+1]
%
%  Reference:
%
%    Louis A. Hageman and David M. Yound
%    Applied Iterative Methods
%    Dover Publications, 2004.
%
%  Inputs:
%
%    real       A(N,N)     -- N-by-N matrix
%    real       xo(N)      -- Initial guess N-vector.
%    real       lambda_ref -- The actual eigenvalue for checking errors
%    integer    maxit      -- The maximum number of iterations.
%    real       tol        -- The error tolerance (defined...?)
%
%  Outputs:
%
%    real       x(N)       -- The computed eigenvector.
%    real       lambda     -- The computed eigenvalue.
%    real       rho        -- The estimated dominance ratio.
%    real       err_norm   -- The norm of the error.
%    integer    iter       -- The number of iterations performed.
%    integer    flag       -- flag (0=okay, 1=maxit reached, -1=bad bad news)
%

if (nargin==6) 
  print = 1;
end
flag    = 1; % unconverged
X_err   = 0;
L_err   = 0;

% Normalize the initial guess.
X  = Xo / norm(Xo, 2);
AX = A * X;
for iter = 1:maxit 
    X     = AX / norm(AX, 2);
    AX    = A * X;
    L     = (X'*AX);  % use Rayleigh quotien estimate for L
    X_err = norm(X - sign(X(1)/X_ref(1))*X_ref, 2);
    L_err = abs(L - L_ref);
    if ( X_err <= tol  && L_err <= tol)
        flag = 0;
        break;
    end;
    if ( mod(iter, 300) == 0 )
        printout(iter,X_err,L,L_err,print)
    end
end
if (flag==0 && print == 1)
    disp(' *** PI: Final Results *** ');
end
if (flag==1 && print==1)
    disp(' *** Warning: PI did not converge ***');
end
printout(iter,X_err,L,L_err,print)

end

function printout(iter,X_err,L,L_err,print)
if (print==1)
  disp(sprintf(...
  ' Iter = %5i, X_err = %6.5e, L = %6.4e, L_err = %6.3e', iter, X_err,L,L_err));
end
end

